library(Seurat)
library(tidyseurat)
library(Seurat)
library(ggpubr)
library(tidyverse)

my.gene <- c('Lepr','Gsdma','Gsdma2', 'Gsdma3','Gsdmc','Gsdmc2','Gsdmc3','Gsdmc4','Gsdmd','Gsdme')

# Bruning 2022 -------------
bruning <- read_rds('hypoMap.rds')

bruning@meta.data |> glimpse()
bruning@meta.data |> count(Author_Class_Curated)
bruning@meta.data |> count(Sex)
bruning@meta.data |> count(Strain)
bruning@meta.data |> count(Author_Condition)
bruning@meta.data |> count(Age)
bruning@meta.data |> count(C25_named)

rownames(bruning) |> str_subset('^Gsdm') |> str_sort()

# major clusters
DotPlot(bruning, features = my.gene, group.by = 'Author_Class_Curated')

# brain spatial sectors
# DotPlot do not allow NA in grouping vars
bruning |>
  mutate(Region_summarized = ifelse(is.na(Region_summarized), 'NA', Region_summarized)) |>
  DotPlot(features = my.gene, group.by = 'Region_summarized') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(x = 'genes', y = 'brain region')

bruning |>
  mutate(Region_summarized = ifelse(is.na(Region_predicted), 'NA', Region_predicted)) |>
  DotPlot(features = my.gene, group.by = 'Region_summarized') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(x = 'genes', y = 'brain region', title = 'region predicted')

DotPlot(bruning, features = my.gene, group.by = 'Author_Region') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(x = 'genes', y = 'brain region')

# diet
DotPlot(bruning, features = my.gene, group.by = 'Diet')

# sex
bruning |>
  filter(Sex %in% c('M','F')) |>
  DotPlot(features = my.gene, group.by = 'Sex') +
  labs(x = 'genes', y = 'Sex')

# age
bruning |>
  filter(!is.na(Age)) |>
  DotPlot(features = my.gene, group.by = 'Age') +
  labs(x = 'genes', y = 'Age')

FeaturePlot(bruning, my.gene)

neur_bruning <- bruning |>
  filter(str_detect(Author_Class_Curated, 'Neuron'))

neur_bruning@meta.data |>
  filter(str_detect(C2_named, 'C2-1')) |>
  count(C25_named)

neur_bruning |>
  filter(str_detect(C2_named, 'C2-1')) |>
  DotPlot(features = my.gene, group.by = 'C25_named')

neur_bruning |>
  filter(str_detect(C25_named, 'C25-8')) |>
  DotPlot(features = my.gene, group.by = 'C286_named') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

# Wu 2021 ----------
wu21 <- data.table::fread('GSE151060_HY_integrated_43261cells_counts.txt.1.gz')

swu21 <- wu21 |>
  column_to_rownames('V1') |>
  CreateSeuratObject()

swu21 |> count(orig.ident)

swu21 <- swu21 |> quick_process_seurat(pcs = 40)

DotPlot(swu21, features = my.gene)

DotPlot(swu21, features = my.gene, group.by = 'orig.ident')

## hpca auto type
hpca <- celldex::HumanPrimaryCellAtlasData()

swu21 <- mark_cell_type_singler(swu21, ref = hpca, new_label = 'hpca_main')

## manual type
marker_wu21 <- c('Vim','Ascl1','Neurog2','Aldh1l1','Slc32a1','Slc17a6','Olig1','Foxj1','Col23a1')
type_wu21 <- c('Radial glial cell','IPC1','IPC2','Astrocyte','GABA-neuron','GLU-neuron','OPC','EC','RGC3')

simpair <- tibble(gene = marker_wu21, wu21_main = type_wu21)

swu21 |> DotPlot(features = marker_wu21)

DimPlot(swu21, cols = DiscretePalette(30), label = TRUE, label.box = TRUE)

swu21 <- FindClusters(swu21, resolution = .5)

res_man <- FindAllMarkers(swu21, features = marker_wu21, only.pos = TRUE) |>
  group_by(cluster) |>
  slice_min(p_val_adj)

res_latent <- res_man |>
  slice_max(avg_log2FC) |>
  left_join(simpair) |>
  group_by(gene) |>
  mutate(wu21_latent = seq_along(wu21_main),
         wu21_latent = str_c(wu21_main, '_', wu21_latent))

swu21 <- res_latent |>
  ungroup() |>
  select(cluster, wu21_main, wu21_latent) |>
  rename(seurat_clusters = cluster) |>
  left_join(swu21, y = _)

DimPlot(swu21, group.by = 'wu21_main', cols = DiscretePalette(10))

DotPlot(swu21, features = my.gene, group.by = 'wu21_main')

# spatial data from www.BrainCellData.org ------
scn.noz <-
  data.table::fread('mission/Duan_GSDMx/Single_Nuc_Cluster_NonZero_Counts.csv.gz')

scn.noz[1:5,1:5]

scn.noz |>
  colnames() |> tail()

scn.gsdme <- scn.noz |>
  select(1, starts_with('Vip='), starts_with('Gsdme')) |>
  as_tibble()

scn.gsdme <- scn.gsdme |>
  mutate(Annotation = str_remove(V1, '.+='),
         gsdme = `Gsdme=ENSMUSG00000029821`,
         vip = `Vip=ENSMUSG00000019772`,
         .keep = 'none')

scn.gsdme |>
  slice_max(gsdme, n = 20)

scn.gsdme |>
  mutate(clus.prefi = str_remove(clus.key, '_.+')) |>
  count(clus.prefi)

scn.meta <- read_delim('mission/Duan_GSDMx/langlieb2023/CellType_Metadata.tsv')

scn.tidy <-
scn.gsdme |>
  left_join(scn.meta) |>
  select(Annotation, gsdme, vip, num_cells_postQC, Max_TopStruct, Max_DeepCCF,
         cell_class)

## allen CCF abbr ------
accf <- read_delim('mission/Duan_GSDMx/allenCCF.csv') |>
  select(name, acronym)

## coarse 12 region --------
scn.coarse <-
scn.tidy |>
  summarise(num.cell = sum(num_cells_postQC),
            gsdme = sum(gsdme) / num.cell,
            vip = sum(vip) / num.cell, .by = Max_TopStruct)

scn.coarse |>
  filter(!is.na(Max_TopStruct) & str_starts(Max_TopStruct, 'NA', negate = T)) |>
  left_join(accf, join_by(Max_TopStruct == acronym)) |>
  mutate(major_region = fct_reorder(name, gsdme)) |>
  pivot_longer(3:4, names_to = 'gene', values_to = 'avg.cnt') |>
  ggplot(aes(avg.cnt, major_region)) +
  geom_col() +
  facet_wrap(~gene, scales = 'free_x') +
  labs(title = 'Gsdme & Vip expression in mouse brain major region',
       x = 'Average expression',
       subtitle = 'Data: Langlieb et al., 2023, Nature') +
  theme_pubr()

## cell types --------
scn.tidy |>
  summarise(num.cell = sum(num_cells_postQC),
            gsdme = sum(gsdme) / num.cell,
            vip = sum(vip) / num.cell, .by = cell_class) |>
  mutate(cell_type = fct_reorder(cell_class, gsdme)) |>
  pivot_longer(3:4, names_to = 'gene', values_to = 'avg.cnt') |>
  ggplot(aes(avg.cnt, cell_type)) +
  geom_col() +
  facet_wrap(~gene, scales = 'free_x') +
  labs(title = 'Gsdme & Vip expression in mouse brain cell types',
       x = 'Average expression',
       subtitle = 'Data: Langlieb et al., 2023, Nature') +
  theme_pubr()

scn.fine <-
  scn.tidy |>
  summarise(gsdme = sum(gsdme), num.cell = sum(num_cells_postQC),
            avg.cnt = gsdme / num.cell, .by = Max_DeepCCF)

scn.fine |>
  slice_max(avg.cnt, n = 20) |>
  filter(!is.na(Max_DeepCCF) & str_starts(Max_DeepCCF, 'NA', negate = T)) |>
  left_join(accf, join_by(Max_DeepCCF == acronym)) |>
  mutate(sub_region = fct_reorder(name, avg.cnt)) |>
  ggplot(aes(avg.cnt, sub_region)) +
  geom_col()
